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ABSTRACT 

Yttria-stabilized zirconia’s high oxygen diffusivity and corresponding high ionic 
conductivity, and its structural stability over a broad range of temperatures, have made the 
material of interest for use in a number of applications, for example, as solid electrolytes in 
fuel cells. At low concentrations, the stabilizing yttria also serves to increase the oxygen 
diffusivity through the presence of corresponding oxygen vacancies, needed to maintain 
charge neutrality. At higher yttria concentration, however, diffusivity is impeded by the 
larger number of relatively high energy migration barriers associated with yttrium cations. 
In addition, there is evidence that oxygen vacancies preferentially occupy nearest-neighbor 
sites around either dopant or Zr cations, further affecting vacancy diffusion. We present the 
results of ab initio calculations that indicate that it is energetically favorable for oxygen 
vacancies to occupy nearest-neighbor sites adjacent to Y ions, and that the presence of 
vacancies near either species of cation lowers the migration barriers. Kinetic Monte Carlo 
results from simulations incorporating this effect are presented and compared with results 
from simulations in which the effect is not present. 

INTRODUCTION 

Zirconia-based materials, and yttria-stabilized zirconia (YSZ) in particular, are of interest for a 
variety of technological applications. Pure zirconia exists in a monoclinic structure below about 
1 100C, with a tetragonal phase stable above that temperature to about 2300C, and a cubic phase 
stable from there to the melting point [1], These phase transitions limit the suitability of zirconia 
for high-temperature applications, particular those that involve frequent thermal cycling. 

The tetragonal and cubic phases can be stabilized via substitutional doping, with aliovalent 
cations such as Y 3+ or Ca 2+ replacing Zr 4+ ions. The cubic phase can be fully stabilized via 
doping with Y 3+ , and the resulting yttria-stabilized zirconia remains in the same cubic fluorite 
phase from room temperature to the melting point. YSZ’s thermal stability and low thermal 
conductivity make it suitable for high-temperature applications such as thermal barrier coatings 
for turbine engine components. 

Oxygen diffusion takes place via diffusive oxygen ion hopping among vacancy sites, with the 
oxygen ion passing between two barrier cations. When zirconia is cation-doped, additional 
compensating oxygen vacancies are formed so as to maintain overall electrical neutrality. The 
addition of these vacancies increases the number of potential hopping sites, increases the 
diffusivity and the diffusive ionic conductivity, which can attain values large enough to make 



YSZ and related materials of interest for use as oxygen sensors, or as solid electrolytes for fuel 
cells. However, the oxygen diffusivity does not increase monotonically with Y 3+ concentration; it 
increases at low concentrations, but reaches a maximum between 8 and 15 mol% Y 2 O 3 , and 
decreases at higher concentrations [ 2 ]. 

In addition, there is evidence for oxygen vacancies to exist preferentially as nearest neighbors to 
either Zr or Y cations. While this issue has remained unresolved, it is potentially important in 
that the presence of an oxygen vacancy in the vicinity of one of the barrier cations may affect the 
value of the migration energy barrier, and thus affect the diffusivity. 

In this work we investigate oxygen diffusion in YSZ using a kinetic Monte Carlo (kMC) 
computer simulation procedure, with barrier energies from density functional calculations. We 
discuss the temperature dependence of the diffusivity, as well as its dependence on yttria 
concentration. We also perfonn density functional calculations of the energetics of the proximity 
of oxygen vacancies to Zr and Y cations, and describe kMC simulations that take into account 
the modification of migration barrier energies due to the presence of these vacancies. 

Fully stabilized YSZ exists in a cubic fluorite structure in which the Zr and Y cations are located 
at sites on a face-centered cubic sublattice, while the oxygen ions are located on a simple cubic 
sublattice whose lattice constant is one-half that of the cation sublattice. Oxygen diffusion takes 
place primarily via the hopping of ions in the [ 100 ] direction to nearest-neighbor vacancies on 
the oxygen sublattice. 

In view of the range of potential applications, considerable experimental and theoretical effort 
has gone into understanding the behavior of YSZ, and, in particular, oxygen diffusion in YSZ 
and similar materials. Computer simulations using a variety of techniques have been performed. 
Schelling et al. [3] have investigated the cubic-to-tetragonal phase transition using molecular 
dynamics simulation, and have correctly predicted the experimentally-observed stabilization due 
to yttrium doping of ZrCfi. Similar work has been carried out by Fabris et al. [4], Fevre et al. 
have investigated the thermal conductivity of YSZ using both Monte Carlo [5] and molecular 
dynamics [ 6 ] techniques. Krishnamurthy et al. have perfonned kinetic Monte Carlo (kMC) 
simulations of oxygen diffusion in YSZ [7] and similar compounds [ 8 ], and have produced 
diffusivities and diffusivity concentration dependence in reasonable agreement with experiment. 
Other molecular dynamics studies of oxygen diffusion in YSZ have been perfonned by by Kahn 
et al. [9], Okazaki et al. [10], Perumal et al. [11] and Shimojo et al. [12]. 

THEORY 

Kinetic Monte Carlo simulation [13-15] is a method designed to study systems where so-called 
“infrequent events” such as diffusive hopping are important. In such systems, molecular 
dynamics simulation, where the accurate representation of atomic vibrations requires a time step 
on the order of femtoseconds, are often inefficient, as a large fraction of computational resources 
are spent simulating the system’s behavior in between the events of interest. By contrast, the 
kinetic Monte Carlo method allows one to concentrate on such events, and to effectively 
consider only the average behavior of the system between such events, while giving up 
infonnation on the detailed atomic trajectories. 



In kinetic Monte Carlo simulations of diffusive hopping, the hopping probabilities are 
detennined by the migration barrier energies, that is, the difference in the total energy when the 
hopping ion is moved from its original position to the saddle point. When the barrier energies are 
known, the hopping rates v AB may be computed from v AB =v° exp (~E AB /k B T) in which v AB and 
Eab are the hopping rate and migration barrier energy for a hop between oxygen sublattice sites 
A and B respectively, and v° is a frequency factor. v° is typically assigned a value between 10 12 
and 10 13 (the value used here) for these materials. For each possible hop, the hopping probability 
Pab can be computed from the hopping rate, with P AB = v AB IT, where T is the sum of hopping 
rates for all possible hops in the computational cell. A catalog of all possible hops, and the 
corresponding hopping rates and probabilities, is created. 

During the kMC process, one of the possible events (that is, a hop defined by the hopping ion, 
the barrier cation species and the target vacancy site) is chosen probabilistically from the catalog 
and executed. Hopping probabilities for all possible hops involving the new vacancy location are 
computed and added to the catalog, while probabilities involving the vacancy's previous location 
are deleted, and the sum of the hopping probabilities is updated. Finally, the simulation clock is 
advanced by a stochastically-chosen time step, with At = -\n(R)!Y, with 0 < R < 1 a random 
number. When the simulation has run for a time t long enough to accumulate statistically useful 
infonnation, the mean square displacement, averaged over all vacancies, is computed. The 
vacancy diffusivity D v is obtained from the Einstein relation (ip 'j = 6 D v t and the ionic 

diffusivity D; is obtained by balancing the number of vacancy and ionic hops, with 
D = C V D V /(I - C v ), where Cv is the concentration of oxygen vacancies. 

The most energetically favorable hopping paths are in the [100] direction, as confirmed by ab 
initio calculations [7]. The hopping oxygen ion proceeds between two barrier cations to a 
vacancy site one-half lattice constant away from the original location. 

RESULTS AND DISCUSSION 

Ab initio energy barriers were computed using the Abinit density functional code[16]. There is 
conflicting evidence that oxygen vacancies preferentially occupy nearest neighbor sites of either 
Y or Zr ions. In order to resolve this ambiguity, we have performed ab initio calculations of the 
total energy of a 2x2x2 cubic YSZ supercell with an oxygen vacancy in a nearest-neighbor 
position with respect to Zr and Y ions. We find that the Y nearest neighbor vacancy position is 
favored over the Zr nearest neighbor position by 0.20 eV. In view of this result, we have. 

Table 1. Oxygen migration barrier energies with and without barrier cation nearest 
neighbor oxygen vacancies. 


Barrier Configuration 

Barrier Energy, eV 

Zr-Zr, no vacancy 

0.706 

Zr-Zr, nn vacancy 

0.568 

Zr-Y, no vacancy 

1.214 

Zr-Y, Zr-nn vacancy 

0.930 

Zr-Y, Y-nn vacancy 

0.925 

Y-Y, no vacancy 

1.941 

Y-Y, Y-nn vacancy 

1.64 




computed barrier energies for all possible barrier ion pairs, with and without nearest neighbor 
vacancies for the barrier cations The results are shown in Table 1. It is evident that there are 
significant differences in the barrier energies with and without a neighboring oxygen vacancy, 
ranging from 0.14 to 0.3 eV. We have therefore included the vacancy-corrected energies in our 
simulations stochastically, based on the relative Zr and Y concentrations. 

The dependence of the oxygen diffusivity on yttria concentration, for a variety of 
temperatures, is shown in Figure 1. The values of the diffusivity are generally consistent 
with experiment [18]. The diffusivity is seen to increase at low yttria concentration, 
reaching a maximum at about 8-25% Y 2 O 3 , depending on the temperature, and decreasing 
at higher concentrations. The concentration at which the maximum occurs is seen to 
increase with increasing temperatures. Both of these points are consistent with 
observation. 
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Figure 1. Oxygen diffusivity versus yttria concentration, incorporating barrier cation 
oxygen vacancy pairing corrections to barrier energies. 

The temperature dependence of the diffusivity is shown in Figure 2, along with results 
from experiment. It can be seen that the values of the kMC diffusivity are reasonably 
consistent with the experimental results of Oishi [18], though the kMC results exhibit a 
smaller range of slopes, and hence a smaller range of activation energies. The kMC results 
generally show an increasing activation energy with yttria concentration, which is 
consistent with the existence of a larger number of high energy barriers containing one or 
two Y 3+ ions. 
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Figure 2. Temperature dependence of diffusivity from kMC results (open icons) and 
experiment (solid icons). 

The concentration dependence of the diffusivity, with and without barrier cation vacancy 
pairing, is shown in Figure 3. It can be seen that in both cases the diffusivity at a given 
temperature is larger for the cases incorporating oxygen vacancies, consistent with the 
lower barrier energies in such cases. The differences are very small for T=1500K, but the 
maximum difference in the Y=2500K diffusivities is approximately fifteen percent. The 
location of the diffusivity maximum with concentration does not appear to change 
significantly when vacancy pairing is included. It is also found that the variation in the 
temperature dependence with and without vacancy pairing is very small. 
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Figure 3. Oxygen diffusivity versus yttria concentration, with (solid icons) and without 
(open icons) barrier cation oxygen vacancy pairing. 





CONCLUSIONS 

We have performed kinetic Monte Carlo computer simulations of oxygen diffusion in yttria 
stabilized zirconia. Isolated oxygen vacancies in a nearest-neighbor configuration with respect to 
a Y 3+ cation are energetically favored by 0.2eV compared with a nearest-neighbor position with 
respect to a Zr 4+ cation. Incorporating oxygen vacancies as nearest neighbors of either of the 
cations in a barrier pair results in a lowering of the barrier energies by 0.14 to 0.3eV, with a 
resulting increase in oxygen diffusivity, though the effect is only significant at high temperatures 
and high yttria concentrations. Oxygen diffusivities computed using barrier energies from 
density functional calculations, with and without vacancy pairing, are consistent in magnitude 
with results from experiment and from other simulations. The dependences of the diffusivity on 
Y concentration and on temperature are also consistent with experiment. 
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